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We propose two ways of estimating the current source density (CSD) from 
measurements of voltage on a Cartesian grid with missing recording points 
using the inverse CSD method. The simplest approach is to substitute local 
averages (LA) in place of missing data. A more elaborate alternative is to esti- 
mate a smaller number of CSD parameters than the actual number of record- 
ings and to take the least-squares fit (LS). We compare the two approaches in 
the three dimensional case on several sets of surrogate and experimental data, 
for varying numbers of missing data points, and discuss their advantages and 
drawbacks. One can construct CSD distributions for which one or the other 
approach is better. However, in general, LA method is to be recommended 
being more stable and more robust to variations in the recorded fields. 



1 Introduction 

A common measure of neural population activity is the local field potential (LFP ), the 
low-frequency part of the extracellular electric potential (iNunez fc Srinivasa 3I2QQ51). The 
LFP is generated by trans-membrane currents in neighboring cells which are usually de- 
scribed on a coarse-grained level by th e current source density ( CSD) ( jPlonseyl Il969l . 
Nicholson &: LlinasI ll97lL iMitzdorj 119851 . iNunez fc Sriniyasan \ mm . In the quasi-static 
approximation the relation of the CSD, C, to the potentials, (j), is 

V(aV0) = -C, (1) 

where a is the electrical conductivity tensor, which for simplicity we assume to be a 
constant scalar (isotropic, homogeneous medium). One consequence of this equation is 
non-locality: is not trivial even in regions where C = 0. This means that the recorded 
LFP may reflect the activity of quite distant cells. 
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When the recordings of LFP at several loc ations are avail able one can attempt re- 
construction of the CSD which generated them (lMitzdorilll985l ). Such recordings can be 
obtained e.g. w ith an electrode with multiple contacts or with a two-d imensional multi- 
electrode array (ICsicsvari et al.ll2003l . iBartho et al.ll2004 lBuzsakill2004l ). 

The simplest method to calcul ate CSD is to use a numerical approxim ation of the 
second derivative of the potential ( Nicholson fc Llinas 1971 . Mitzdori llQSSi ). e.g. in case 
of 1-D electrode with equidistant contact points spaced by h one obtains (for interior 
contacts): 

. . (j){zi + h)- 2(f){zi) + (j){zi - h) 
C-iZi) = -a — . (2) 



Such an approach has several disadvantages. One of them is that Eq. (I2|) cannot be applied 
to boundary points. This is particularly inconvenient in case of t wo- or three-dime nsional 
data, where the boundary may comprise majority of the points (IL^ski et al.ll2007ll. 



A nother method for estimating the CSD is the inverse CSD (iCSD) method (iPettersen et al 

2006 . L^ski et aD 200?! ). Here one does not try to use Eq. directly (which is the case 



in traditional CSD). Instead, the idea is to establish a one-to-one relation F between 
measured voltages and CSD distributions via inversion of the forward solution. This is 
achieved in the following way: assume N recording points on a Cartesian grid (one-, two- 
or three-dimensional). Consider A^-parameter family of CSD distributions — this means 
that given the values of the N parameters one can assign a value of CSD to each spatial 
position. Then the values of the potential, 0, on the grid can be obtained by solving 
a well-posed boundary value problem related to the elliptic partial differential equation, 
Eq. ^ (forward solution). Therefore, the N measured voltages are functions of the CSD 
parameters. If the family and the parameterization are chosen well, one can invert this re- 
lation and from the N measured potentials recover the N parameters of CSD. Usually one 
parameterizes the CSD with its values on the measurement grid and inter polates between 



the grid points, linearl y or with splines, but there are more possibilities (IPettersen et al. 



2nnalL(;ski et al.ll2nn7h . 



In the reconstruction one may assume that the CSD is non-zero only inside the grid. 
Usually, however, the actual CSD extends in the tissue beyond the grid set by t he measure 



ment points and such an assumption would lead to large reconstruction errors (ILeski et al 



20071 ). To avoid them one may use a trick of extending the grid spanning the CSD by 
an additional layer. One can set the CSD values at the additional no des to zero or du - 
plicate the value from the neighboring node of the original grid. In (IL^ski et al.l 120071 ) 
these two approaches were denoted by B or D boundary conditions, respectively, and it 
was shown that they improve the reconstruction quality. Note that the new CSD family 
is still parameterized with its N values at the original grid and that the intention of such 
a procedure is to improve the reconstruction fidelity inside the grid and not to estimate 
the CSD outside the grid. 



2 Inverse CSD on incomplete data 

A practical problem in the application of iCSD method to real datasets is how to deal with 
missing recording points. Such cases arise surprisingly often in real experiments. There 
may be several reasons for this: a contact of a multi-electrode may not be functioning, 
some channels may be used for other purposes (e.g. stimulation), or the experiment may 
be terminated early before all the data are collected. 
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One way to deal with such data is to patch them with the mean of the neighboring 
potentials of a missing contact (IL^ski et al.ll2008l ) . We denote this method by LA for local 
averages. This approach means replacing the missing true potential at a point with a 
linear approximation estimated from the neighbors. 

A more elaborate alternative is to reduce the size of the CSD grid and find the least- 
squares solution to such an overdetermined system. That means choosing such values of 
the parameters of the CSD spanned on a smaller grid which minimize the sum of squared 
differences in potential at all the available electrode points. We denote this method by LS 
for least squares. The advantage of this approach is that we use only the available data 
without making any assumptions about the missing recordings, so it seems to be better 
motivated than LA. However, we decrease the spatial resolution of the reconstructed CSD 
so it is hard to tell a priori which method is better. 



3 Gaussian sources 

To find out which of the proposed approaches works better we first tested them both on 
three-dimensional surrogate Gaussian sources (Fig. [1]). We calculated the potentials on a 
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Figure 1: A) Gaussian sources studied in Section [3l four consecutive slices (x = 1 . . .4) 
through the volume. Electrode positions are marked with circles. B) Reconstruction of the 
CSD distribution spanned on the full 4 x 10 x 4 grid from the set of potentials calculated 
at the nodes of the grid (denoted by x's). C) Reconstruction using LS method on complete 
data; spanned on a smaller, 4x8x4 grid. 
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grid of (x, 2) G 4 X 10 X 4 equally spaced points. Details of the structure of the sources 
and calculation of LFP are given in the Appendix. This choice of the sources and the grid 
was motivated by a recent experimental study of evoked potentials in the barrel cortex of 
the rai0. The sources were elongated along the y axes so that the conclusions would hold 
for cortical dipoles generated by active pyramidal cells. 

For the tests we removed a number of virtual 'recording points', reconstructed the 
CSD u sing both LA and LS methods, and compared the normalized reconstruction 



errors (ILeski et al.l 120071 ): e = f{C — C)^dx/ f C^dx, where C is the original and C is 
the reconstructed CSD. For the LS method we used a grid of 4 x 8 x 4 points which 
covered the whole space occupied by the original grid. This implied larger spacing in 
the y direct i on. T he iCSD reconstruction was performed with the Matlab scripts from 
L^ski et al.l ( 120071 ) . modified for the situation at hand. We used not-a-knot splines with D 
boundary conditions!^ 
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Figure 2: Comparison of local averages (LA) and least squares (LS) methods of recon- 
structing CSD from incomplete data. A) Results of reconstruction in the case of a single 
recording point removed from the grid. Normalized reconstruction error for all 160 pos- 
sibilities, LS: x's, LA: o's, sorted according to the LS error. B) Histogram of normalized 
reconstruction errors for a pair of grid points removed. Thick bars: LA, thin bars: LS. 
Inset: outliers in the LS method. C) Comparison of LA (o's) and LS (x's) methods for 
varying number n of recording points removed from the grid (X axis). Y axis: average 
logarithm of normalized reconstruction error, error bars are ± standard deviation, for the 
best 90% out of 2000 random choices of removed points (except n = 1 where 90% of all 
160 possibilities are taken). D) Same as C, but for the worst 10% of the cases. 

Standard iCSD reconstruction from the calculated potentials gives reconstruction error 



^ J. Ka miriski, private com munication. Note that there are more sophisticated models of the cortical 
LFP, e.g. i Tenke et al. 1993), but for these tests we found our simple model adequate. 



^By D boundary conditions we mean solution on a larger grid with one extra layer added in every 
direction beyond the original. We assume identical CSD at the added layer and its nearest neighbor in the 
original grid. 'Not-a-knot' splines are the cubic splin es implemented in Matlab, they diSer from 'normal' 



spHnes in the conditions at the extreme points. See lL^ski et al.l l|2Q07l ) for details 
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e of 0.14%. This indicates the quality of the reconstructed approximation of the smooth 
Gaussian sources by a set of splines spanned on the grid of recorded points. The LS 
method applied for the reduced grid gives e = 0.21% which shows how little information 
is lost when the number of nodes of the grid used for the reconstruction is reduced by 
20%, see Figure [TJ This is possible in this case thanks to the relatively large extent and 
slow variation of the sources in the y direction: the sparser grid is still dense enough to 
effectively sample the sources. In general, the degradation of the reconstruction quality 
caused by using a sparser grid will strongly depend on how rapidly the CSD varies in 
space, see Section [H 

We have scanned all the 160 cases of one recording point withdrawn. The LS method 
gives stable reconstruction error from e = 0.21% to e = 0.26%. The error of the LA 
method ranges from e = 0.14% (which means that the missing datum was indeed the 
mean of its neighbors) up to 2.1%, see Fig. [2]A. 

There are 12720 possible choices of a pair of electrodes for the set of recording points 
considered. Here the results are more intricate: as before, LS typically gives smaller errors 
(Fig. [2]B), but from time to time a huge error occurs, with e reaching 270% (63 outliers 
shown in an inset in Fig. [2]B). The outliers can occur only for specific configuration of 
the missing pair (xi, Zi), (x2, ?/2, -22)- The necessary (but not sufficient) conditions for 
large errors are (xi, zi) = (xa, ^2) and (1/1,1/2) E {(1, 2), (1, 3), (2, 3), (9, 10), (8, 10), (8, 9)}. 
There are only 96 of such troublesome pairs and all the outliers in Fig. [2j3 are of this type. 

With growing number of missing recording points the reconstructions become less and 
less reliable with the LS method becoming monotonically worse with respect to the LA 
method. Fig. [2]CI, D. Interestingly, the distributions of the results qualitatively have the 
same character as in the case of two missing electrodes: that is, the errors of LA method 
have a unimodal distribution while the distributions of errors in the LS approach have two 
modes, one with results better than for LA, the other with extremely large errors. The 
mean error of the LA method also grows but huge errors do not occur. Figures ^P, D 
show the results obtained for 2000 random choices of the missing recording points. 



4 Experimental data 

The second part of the test of the two methods was performed on three-dimensional record- 



mgs 



in the rat forebrain of potentials evoked by the defiection of a bunch of whiskers (iLeski et al 



20071 ). The recordings were made on a grid of 4 x 5 x 7 points. Here we analyse t he same 



two representative latencies which were used as illustrations in iL^ski et al.l (120071 ) , where 
the dataset is described in detail. 

We perform the same analysis as for the Gaussian sources, apart from the fact that now 
we do not know the real CSD. Therefore, as the reference C we take the reconstruction 
spanned on the full 4x5x7 grid calculated from the full set of recordings. Such C is the 
best representation of real CSD in the tissue available to us. 

Figure [HK shows the reference data set and Fig. [3|3 shows the CSD reconstructed by 
LS method on a sparser, 4x5x6 grid, from the complete set of recorded potentials, 
3.5ms after the stimulus onset. Already here we can observe how the intricate structure of 
activation in the tissue is distorted (e = 21%) when using a smaller spanning grid which 
was not the case for the Gaussian sources modeling the cortical CSD (Sec. [3l). 

Clearly, performing reconstructions from incomplete data, we expect the distortions to 
grow. Results of the test of the two methods are shown in Fig.Sl Fig. [HA shows the results 
of the iCSD reconstructions from data with one electrode removed. As in the previous 
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Figure 3: Reconstructions of CSD from experimental data, t = 3.5ms after stimulation of 
vibrissa. Each row presents a three-dimensional region of the rat forebrain. The electrode 
positions (4x5x7 grid) are marked with circles, nodes of the grid are marked with 
x's. A) The reference data set: CSD reconstructed on the full electrode grid. B) CSD 
reconstructed from the complete set of recordings but spanned on a sparser 4x5x6 grid. 
Note that some sources are not adequately sampled using the sparser grid. 
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case, the distribution of errors of the LS method is bimodal with very narrow modes, while 
the distribution of errors of the LA method is rather broad (Fig. [4|3)- However, unUke the 
previous case, the LA method is almost always better. This difference is preserved as the 
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Figure 4: Comparison of LA and LS methods of reconstructing CSD from incomplete data, 
see caption of Figure El The data used here are the same as in Figure [3l 

number of removed points is increased. For both methods the mean error of reconstruction 
grows with the number of removed points, which is expected. However, the distribution 
of errors for the LS methods gets wider, while the distribution of errors for LA method 
gets more narrow. This is true for both the best 90% cases and for the 10% worst. In 
practice this means that the LS method for such complicated CSD distributions is not 
recommended. 

Such behavior was typical for this data set for the time frames we inspected. For 
illustration and comparison we show the reconstructions from the complete data on the 
original (Fig. EK) and smaller (Fig. [5j3) grids as well as the results of the same analysis 
for the recordings taken 15ms after the stimulus onset (Fig. [6]). 



5 Discussion 



Reconstruction of the current source density generating recorded extracellular potentials 
from these potentials is an ill-defined problem. The reason is that there is an infinite 
number of different distributions which could lead to the same recordings. Nevertheless, 
as the CSD is much more local reflection of the neural activity than the potentials, there 
were r nanv attempts to find a viable reconstruction of the sources from the measured 
fields (INicholson &: Llinaslll97ll . lMitzdorilll985l . IPettersen et al.ll2006l ). One recent candi- 
date whi ch has a number of advantages oyer th e classical approach is the inverse CSD 
method (IPettersen et al.l l2006l . iL^ski et al.l 120071 ). It has been originally developed for 
situations where a set of recordings was collected on a regular rectangular grid. Given 
the construction of the method it is unclear how to proceed when one of the recordings 
is missing, due to a failure of one of the electrode contacts or in other cases. We have 
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Figure 5: Reconstructions of CSD from experimental data, t = 15ms after stimulation of 
vibrissa. For description see caption of Fig. [3l 
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Figure 6: Comparison of LA and LS methods of reconstructing CSD from incomplete data, 
see caption of Figure [2l The data used here are the same as in Figure [5l 
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discussed here two approaches which might enable the apphcation of iCSD method to sets 
with incomplete data. 

Local Averages method (LA) is simple, stable, and the results are never very bad 
(normalized error of the order of a few percent), even for a relatively large number of 
missing data points. The Least Squares method (LS) seems attractive as it does not 
assume anything about the missing data. The distribution of errors is usually bimodal 
with two narrow modes. Usually, the errors are within a small range dominated by the 
effect of the sparser grid, however, for a subset of cases which is growing with the number 
of missing data, the errors can be extremely large. 

The respective quality of reconstruction for the two methods depends on the structure 
of the original sources, and (especially for LS) on the specific location of the removed 
points. Our tests on the sources modeling the dipole distributions of the cortex (Section [3]) 
with the grid shrinking along the dipole show that for a small number of missing recording 
points (<5) the LS method usually gives smaller errors than the LA. However, for more 
complex thalamic sources (Section S]) the LA method is usually far better for any number 
of removed points. 

A priori it is not obvious which method to choose in analysis of experimental data, 
when the original CSD is unknown and is to be found. Our recommendation is to use the 
LA method in all cases. Despite its simplicity it seems to be more stable and leads to 
smaller errors, especially for complex distributions, thus it becomes our method of choice. 
If the potential seems to vary relatively slowly along one direction of the grid and the 
missing data are not nearest neighbors lying at the edge, the LS method might also be 
worth trying, but in general we do not recommend it. 

One may wonder if it is possible to improve the technique beyond the proposed ap- 
proaches. One way would be to consider CSD distributions spanned on the available 
recording points which would not necessarily form a full regular grid. However, this seems 
rather difficult to implement in full generality, as the spline coefficients would have to be 
calculated from the scratch for every distribution of the recording points, and the matrix 
connecting the potentials with the CSD parameters would have to be calculated for ev- 
ery distribution adding substantially to the computational overhead. A more promising 
approach seems through the application of statistical methods. For example, one way we 
plan to follow in the future is to use an overdetermined basis of Gaussian sources and 
search for efficient projections of the recordings on this basis. 

Appendix: Gaussian test sources 

The Gaussian sources used in the test in Section [3] were of the form 

jx - Xj^ _ (y - _ {z- Zjf 
2{aTY 2{atY 2{aTY _ ' 

with the coefficients given in Table [H Fig. [T] shows four parallel sections of the sources 
which together pass through all the nodes of the grid (in the region spanned by the 
virtual recording grid). To calculate the potentials we truncated the sources to the region 

-1 < a; < 6, -1 < y < 12, -1 < z < 6. 
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Table 1: Coefficients of the Gaussian sources. Origin of the grid is {x,y,z) = (1, 1, 1). 
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